# Function to generate plots for the paper
genPlot <- function(mydata, DV, dvName, ylim, t1, t2, run_var = 'run_var'){
	mayData <- mydata

	m1 <- rdrobust(y = mydata[, DV]
			, x = mydata[, run_var]
			, p = 1
			, c = 0		
			, level = 95 
			, all = TRUE)
	p1 <- rdplot(y = mydata[, DV]
			, x = mydata[, run_var]
			)

	par(mar = c(5, 5, 2, 2))
	plot(x = p1$vars_bins$rdplot_mean_x[p1$vars_bins$rdplot_mean_x > -55 & p1$vars_bins$rdplot_mean_x < 55],
		y = p1$vars_bins$rdplot_mean_y[p1$vars_bins$rdplot_mean_x > -55 & p1$vars_bins$rdplot_mean_x < 55],
		ylim = ylim,
		xlim = c(-50, 50), 
		pch = 20, 
		col = 'grey',
		xlab = 'Difference in Vote Share, Municipal Election (t-1)',
		ylab = paste0(dvName, ' (t)'),
		axes = FALSE,
		cex.lab = 2
		)
	lines(x = p1$vars_poly$rdplot_x[p1$vars_poly$rdplot_x < 0]
		, y = p1$vars_poly$rdplot_y[p1$vars_poly$rdplot_x < 0], lwd = 2)
	lines(x = p1$vars_poly$rdplot_x[p1$vars_poly$rdplot_x > 0]
		, y = p1$vars_poly$rdplot_y[p1$vars_poly$rdplot_x > 0], lwd = 2)
	axis(1, at = c(seq(-50, 50, by = 25)), cex.axis = 2)
	axis(2, cex.axis = 2, at = round(seq(ylim[1], ylim[2], length.out = 5), 1))
	abline(v = 0, lty = 2)
	tauvalue = round(m1$Estimate[1, 'tau.bc'], 3)
	text(x = -30, y = t1, paste0('Effect = ', tauvalue), cex = 2.2)
	text(x = -30, y = t2, paste0('p-value = ', round(m1$pv[3, 1], 3)), cex = 2.2)
}



#----------Sensitivity Test ----------#

# Function to run sensitivity
sensitivity <- function(DV, ylb, yub, mydata, run_var, ylabel){
  betas <- c()
  lb <- list()
  ub <- list()
  h <- seq(0.5, 100, by = 0.05)
  for(i in 1:length(h)){
    out <- rdrobust(y = mydata[, DV]
                    , x = mydata[,run_var]
                    , h = h[i]			
                    , c = 0
                    , level = 95
                    , all = TRUE)	
    betas[i] <- out$Estimate[,'tau.bc']
    lb[i] <- out$ci[3, 1]
    ub[i] <- out$ci[3, 2]
  }
  
  # Create Plot
  par(mar = c(5, 7, 2, 2))
  plot(x = h, y = betas, type = 'l', 
       ylim = c(ylb, yub), xlab = 'Bandwidth',
       ylab = ylabel,
       cex.lab = 1.75, axes = FALSE)
  lines(x = h, y = unlist(lb), lty = 3, lwd = 1.25)
  lines(x = h, y = unlist(ub), lty = 3, lwd = 1.25)
  abline(h = 0, lty = 2, lwd = 0.5)
  axis(1, cex.axis = 1.5)
  axis(2, cex.axis = 1.5)
}

#----------Difference in Coefficients ----------#

### This function is used in the replication to run a Z-Test
zDif <- function(beta1, beta2, se1, se2){
  dif <- beta1 - beta2
  se <- sqrt(se1^2 + se2^2)
  zscores <- dif/se
  pvalues <- ifelse(zscores > 0, (1 - pnorm(zscores)) * 2, pnorm(zscores) * 2)
  out <- data.frame(difference = dif, SE = se, zscore = zscores, pvalues = pvalues)
  return(out)
}
